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The elimination procedure as described by Williams has 
been coded in LISP and FORMAC and used in solving systems 
of polynomial equations. It is found that the method is very 
effective in the case of small systems, where it yields all solu¬ 
tions without the need for initial estimates. The method, by 
itself, appears inappropriate, however, in the solution of large 
systems of equations due to the explosive growth in the inter¬ 
mediate equations and the hazards which arise when the 
coefficients are truncated. A comparison is made with difficulties 
found in other problems in non-numerical mathematics such as 
symbolic integration and simplification. 

Let us oonsider the following problem: 

Given a system of polynomial equations with real 
coefficients, 

Piixi , *2, • ■ • , x„) = 0, i = 1,2, ■■■ ,n, 
find the solutions (common zeros) of the system (assuming 
the consistency of the system). 

The standard attach of numerical analysis on this 
problem is by an iterative method based on Newton’s 
method or the method of steepest descent. These methods 
assume a good estimate for the solutions and they may, 
in fact, fail to converge otherwise. Such an estimate will, 
in general, be required for each desired solution. Obtaining 
the estimates can be a very time-consuming process and 
obtaining sufficiently good estimates may be a prohib¬ 
itively costly operation. 

The approach considered here is the one taken by the 
algebraists of the last century which led to the elimination 
method. The elimination method reduces the problem to 
that of finding the zeros of a polynomial in a single variable. 
It is very similar to the Gaussian elimination method used 
in solving linear equations. 

In the Gaussian elimination procedure (Figure 1) the 
variable Xi is eliminated in the step labelled A. The 
corresponding step B for the polynomial case (Figure 2) is 
a call to a subroutine (Figure 3) which produces two 
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poljmomials which we shall call the divisor and the elim- 
inant. The eliminant is independent of Xi and the divisor 
is usually of lower degree in x, than the original equations. 
Note the .similarity of the form of the transformation 
involved in step C of the subroutine and that in step A. 


for i ■.= 1 stop 1 until Ji—1 do 
for J := t -t- 1 step 1 until n do 
A Pi •■= aaPj — ajiPi 

Fig. 1. Gaussian elimination for linear equations 

Pi = ^ atiXk + oe.M.1 = 0, Pi = YL 0 


for t := 1 step 1 until n-1 do 
for j := f-Fl step 1 until n do 


Pi := divisor of Pi and Pi with respect to an 
Pi eliminant of Pi and Pi with respect to an 


Fig. 2. Elimination for polynomial equations 


C |—>0 = S.n - Rr&xi • 

q = degree of Q in a:i 
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While step A completely eliminates the variable, step C 
can only eliminate the highest power of the variable 
appearing in the two equations. In the linear case the 
transformation in A is best done by an algebraic processor 
such as Foktran. However, in the nonlinear case a 
processor with polynomial manipulation capabilities is 
essential in order to efficiently perforin the transformation 
indicated in step C. We shall call the resultii^; system of 
equations, P,-, i = 1, 2, • • • , n, the reduced system. 

The following theorem [1] provides the justification for 
the elimination procedure used in the nonlinear case. 

Theorem. A common zero of Pi and Pj is a common zero 
of the eliminant and divisor of Pi and Pj loUh respect to x, . 

Hence, in order to find the solutions to the original 
system of equations wc need find only the solutions to the 
reduced system. In the reduced system each equation P,- 
is independent of Xj if j < i. Typically the reduced system 
also has the following character: The reduced equations 
Pi ,i = 1, 2, • ■ • , a —1 are linear in variable Xi , and 
equation P„ is of high degree in Xn , where the degree 
reflects the number of solutions of the system. Thus P„ 
can be solved (in principle) by a root-finding program and 
the complete solutions can be obtained by successively 
solving P„-i, • • • , Pa, Pi for a:„-i, • • • , *2, aii as equa¬ 
tions in one unknown in a backward substitution process 
similar to that used in the Gaussian elimination procedure. 

Although the above is the usual case there are some 
exceptions. For one, the first n—1 equations of the reduced 
system need not be linear. This complicates the backward 
substitution somewhat, but normally poses no real prob¬ 
lem. Another possibility is that solutions of the reduced 
system may be extraneous. This is usually noted in the 
backward substitution phase, when an equation of the 
reduced system becomes identically isero. Such solutions 
are then ignored. Finally, we may find that one of the 
eliminants generated during the elimination process will 
vanish identically. This occurs when equation P,- and Pj 
have a factor in common involving Xi. This difficulty may 
sometimes be overcome by rearranging the equations and 
starting the process over again. Otherwise, the equations 
can bn divided by the common factor, thereby generating 
subsystems of equations each of which must be solved. This 
difficulty is troublesome, but fortunately it docs not occur 
often in practice. 

Williams [1] described how the complete procedure could 
be programmed for a digital computer using a polynomial 
manipulation system, but he lacked access to a sufficiently 
large machine on which to cany out his ideas. The pro¬ 
cedure has been programmed both in Lisp [2] and Foemac 

[3] for the IBM 7094. An example of the solution of a 
small system is given in Figure 4. The solution of another 
system is described in Figure 5. The program uses single 
precision arithmetic throughout except when finding the 
roots of equations. The latter process is performed in 
double precision with a program using Muller’s method 

[4] . It is interesting to note that Bareiss’ root-finding 


program [a] makes use of elimination in a limited context 
in extending the Gracffe root-squaring technique. 

It a system has more variables than equations, ive call 
the variables which are not eliminated parameters of the 
system. J. Israel has written a program in F oemac which 
obtains solutions in closed form to systems containing 
parameters in which the equations in the reduced system 
arc of degree 4 or less in their respective variables. The 
solutions (in terms of the parameters) arc obtained from 
the formulas which are available in such cases. 

The experience gained with the programs leads to the 
following conclusions. 

(1) The method yields fast and accurate results for 
small systems of equations (e.g., three quadratics). No 
starting values are needed and all the solutions (complex, 
in general) are obtained except in degenerate ca.ses. 

(2) When the reduced system is of degree 4 or less in 
each of its variables, standard formulas can bo used to 
obtain solutions in closed form for the original system in 
terms of its parameters. 

(3) For slightly larger systems (e.g., three quarries) 
the degree of the final equation in the reduced system may 
be so large that one cannot obtain sufficiently accurate 
solutions for it with existing root-finding processes. With 
an estimate for a zero of the reduced equation F» such a 
program can find any single desired solution. Backward 
substitution yields the complete solution. In such a case 
the elimination method is superior to existing numerical 
procedures in that it requires an estimate for only one 
variable instead of all. 


C = 1-A-l-B x-Zv^ + v + l = ti 

D = 1-B-l-C-a: (.Zy* - y -‘Ifx+ 2y‘- y - \ = 

E = (i>p-y-2)G -1-D gy - Zy‘ - + 2y + I = Z 

Reduced eyeteni 


Solutions o/E: y = 1, —4, Vi, — Vi 

SuhsttLutton m Cr x = 1, —i, —Vi, Vi 

Time for solution with Formac program: approx. 8 se 


Reditced system 
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(4) For large systems of equations, the final poly¬ 
nomial in the reduced system can be of such high degree 
that it will be virtually impossible to solve for all the roots, 
and the solution for a single root will face a roundoff error 
problem. In many situations, however, the large final 
polynomial can be avoided by overdetermining the system, 
i.e., by adding more equations without any additional 
unknowns. A nontrivially overdetei'mmed system usually 
has far fewer solutions than the original system. The price 
for obtaining a final equation of small degree is the added 
size of the intermediate polynomials and the possible 
increased error due to roundoff. 

We find that there are two types of limitations to the 
computer implementation of solutions to problems in 
non-numerical mathematics. One is due to recursive un¬ 
solvability, the other' is due to the explosive growth in the 
intermediate variables in recursive problems. It is a 
familiar fact that deciding the theorems in tlie lower 
predicate calculus is a recursively unsolvable process and 
hence necessitates the use of heuristics in tlie general 
case. Biohardson [6] has shown that tlie problem of 
matching elementary functions [7] is also recursively un¬ 
solvable. The elementary functions contain the usual 
arithmetical, trigonometric, exponential and logarithmic 
functions of real variables, and the matching problem asks 
whether an elementary function is equivalent to zero. 

It is interesting to note that -the difficulties found in 
matching programs such as Martin’s [8], and Tobey’s [3] 
occur with the introduction of the same functions (log and 
square root) that are critical in Kichardson’s proof. 

Simplification has properties which make it, in general, 
an even more insidious problem than matching. The 
criteria for simplicity are not fixed but may vary with the 
situation. For instance, it may be easier to integrate 
4a:®[a:V((^*)®+l)l than the eqnivalnnT, and apparently 
simpler 4a:V(a:‘^-fl). Suppose that the criteria could be 
fixed (e.g., zero is the simplest form of a function equiv¬ 
alent to it), then simplification behaves like a minimiza¬ 
tion problem in a iioiie-too-well structured space. It is 
difficult to tell whether one has the minimum (i.e., sim¬ 
plest) expression, or what transformation to apply to get 
do'ivii to the minimum. In fact, the inininiuni caimot alwa.ys 
be recognized in a finite number of steps because by 
finding it we could solve the matching problem. 

Bichardson’s unsolvability result for matching also 
carries over to the integration problem for the elementary 
functions. It is hard to say how much Slagle’s integration 
prograni [7] suffers from this difficulty, but the result 
certainly explains the difficulties mathematicians have 
had in solving problems in integration and differential 
equations. 

The limitations of procedures designed to perform 
manipulations (e.g., inversion) on matrices wdth non¬ 
constant entries are ccrtainlj’- not due to recursive un¬ 
solvability. The difficulties in these problems are due to the 
highly nonlinear growth of the intermediate and final 


results, in general. Similarly, the difficulties wo are faced 
with in the solution of systems of polynomial equations 
are due to the explosive gron’th in the degree of inter¬ 
mediate equations and their coefficients in large systems. 

Further remarks about the difficulties due to the 
coefficients are in order. If the coefficients of the original 
system of equations are not rational (hence introducing a 
starting error in the computer realization of the solution), 
then the instability of the problem will lead to gross errors 
in the solution of large systems. Suppose we truncated the 
coefficients in some stage of the solution in order to gain 
space and time, then the solution for sufficiently large 
systems will be in great error due to the roundoff error 
introduced coupled with the instability of the problem. 
Hence in order to obtain meaningful solutions in the 
elimination of large systems, the coefficients should be 
exact at all stages. In Williams’ paper these difficulties 
seem to have been underestimated. 

Recently Collins [8] obtained some major results in this 
area. He derived a small bound on the growth of the length 
of the coefficients in the problem of finding the greatest 
common divisor (GCD) of two polynomials. Usually GCD 
algorithms are essentially extensions of the Euclidean 
algorithm to the polynomial case, and are similar to the 
subroutine m Figure 3. Collins shows that the gro^vth of 
the coefficients is linear in the degree of the original equa¬ 
tions and the length of their coefficients, and is not nearly 
as large as was previously estimated. He furthermore has 
programmed a method for obtaining a reduced system to 
a system of polynomial equations which damps the growth 
of the decree of the intermediate equations and their co¬ 
efficients. We believe these results constitute a break¬ 
through in the solution of this problem by the techniques 
of algebraic manipulation. However, we doubt that any 
single procedure can change the exponential character of 
the problem. One should keep in mind that the number of 
solutions and their absolute values are, in general, exponen¬ 
tial functions of the coefficients, degrees and number of 
equations of the original system. While ideally one would 
like the degrees of each of the polynomials of the reduced 
system to be as small as possible (the danger then being 
the accumulation of errors during the backward substitu¬ 
tion, forcing the solution of polynomial equations with 
inexact coefficients), it seems unlikely that one can avoid, 
in general, simulating a large part of the search process 
necessary to locate the solutions in an a-dimeiisional 
complex space. In numerical analysis this search is 
explicit; the impheit methods in algebraic manipulation 
tend to lead to increasingly large intermediate variables.' 
We venture a guess that Kutta’s system [10] of 16 equa¬ 
tions in 15 unknowns (used in constructing a fifth-order 

' These are systems of equations whose reduced system must 
contain exponentially large polynomials. Consider the system 

contain a polynomial divisible by one of the polynomials — 
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RungeJvutta method), which has, to the best of our 
knowledge, withstood solution by the techniques of 
numerical analysis, is going to remain unsolved by the 
techniques of algebraic manipulation. 

Although there exists for the above two types of prob¬ 
lems in algebraic manipulation no single program which 
will handle all po.saiblc inputs efficiently or at all, there 
are many special categories where such programs can and 
do exist. Even though simplification in general is not 
possible, the simplification of rational functions reduces to 
the GCD problem which is recursive. The integration of 
rational functions in closed form is performed by a program 
written b3^ Manovc [11]. Surprisingly, the integration 
problem for the algebraic functions (these include rational 
function.s and roots of rational functions) appears to be 
recursive also. Hardy [12] conjectured that this problem 
could not be solved in a finite number of steps. 

For the recursive problems which experience growth of 
the intermediate variables there exist similar categories. 
One category is composed of those histances which are so 
simple that reasonable algorithms do not manage to blow 
them up. The simplicity of the algoritlim to be used in 
solving more difficult cases need not be a virtue. For 
example, in elimination one profits from factoring the 
equations because this inhibits their growth. This is the 
manner in which a high school student would proceed on 
the system in Figure 4. The Alpak algorithm for the 
GCD problem [13] makes use of special oases of factoring. 
A more general factoring algorithm is implemented in 
Engloman’s system [11]. Modifications can also be made 
to the basic elimination algorithm to prevent the introduc¬ 
tion of any extraneous roots. In this connection it should 
be noted that Williams did not account for all the extra¬ 
neous roots introduced by the elimination procedure. 

Finally, each particular instance of a problem can be 
considered as a special oategorj’- since the user may have 
some insight into it. Symbolic manipulation systems 
should have sufficient flexibility for the user to translate 
his insight into an algorithm simply. .Since a user fre¬ 


quently cannot foresee the possible consequences of a 
proposed transformation, a system with capabilities for 
man/maehine interaction has many advantages. Martin’s 
displaj"^ program for mathematieal expressions [14] is a 
large step in this direction. We believe that difficult 
problems can be most efficiently solved, if they can be 
solved at all, on systems possessing such capabilities. 
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Informal Evening Sessions 

One of the features of the Sjunposium was a procedure copied from Share, in which in¬ 
formal evening discussions were suggested and led by attendees themselves. Sessions were 
held on the following topics: 

Implementation Techniques 

Approximations to Solutions of Differential Equations 

Mathematics and Metamathematics of Algebraic Manipulation Languages 

LISP II 

Natural Languages Processing 

Symbol Manipulation Techniques As Applied to Graphical Languages and Topological 
Models 

FORMAC 360 External Specifications 

L« — A Computer-Made Movie Showing a Programming Example in L' 

It is worth noting that the LISP II meeting was attended by over 50 people, and the 
session on Metamathcmatics by a vigorous group of about 65 people, even though the time 
was 6-8 p.m.—J. E. S. 
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